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Classical age-structured mass-action models such as the McKendrick-von Foerster equation have 
been extensively studied but they are structurally unable to describe stochastic fluctuations or 
population-size-dependent birth and death rates. Stochastic theories that treat semi-Markov age- 
dependent processes using e.g., the Bellman-Harris equation, do not resolve a population’s age- 
structure and are unable to quantify population-size dependencies. Conversely, current theories 
that include size-dependent population dynamics (e.g., mathematical models that include carrying 
capacity such as the Logistic equation) cannot be easily extended to take into account age-dependent 
birth and death rates. In this paper, we present a systematic derivation of a new fully stochastic 
kinetic theory for interacting age-structured populations. By defining multiparticle probability 
density functions, we derive a hierarchy of kinetic equations for the stochastic evolution of an 
ageing population undergoing birth and death. We show that the fully stochastic age-dependent 
birth-death process precludes factorization of the corresponding probability densities, which then 
must be solved by using a BBGKY-like hierarchy. However, explicit solutions are derived in two 
simple limits and compared with their corresponding mean-field results. Our results generalize both 
deterministic models and existing master equation approaches by providing an intuitive and efficient 
way to simultaneously model age- and population-dependent stochastic dynamics applicable to the 
study of demography, stem cell dynamics, and disease evolution. 


INTRODUCTION 

Age is an important controlling feature in populations of living organisms. Processes such as birth, death, and 
mutation are typically highly dependent upon an organism’s chronological age. Age-dependent population dynamics, 
where birth and death probabilities depend on an organism’s age, arise across diverse research areas such as demog¬ 
raphy [1], biofilm formation [2], and stem cell proliferation and differentiation [3, 4]. In this latter application, not 
only does a the cell cycle give rise to age-dependent processes [5, 6], but the often small number of cells requires 
a stochastic interpretation of the population. Despite the importance of age structure (such as that arising in the 
study of cell cycles [5-7]), there exists no theoretical method to fully quantify the stochastic dynamics of aging and 
population-dependent processes. 

Past work on age-structured populations has focussed on deterministic models through the analysis of the so-called 
McKendrick-von Foerster equation, first studied by McKendrick [8, 9] and subsequently von Foerster [10], Curtin and 
MacCamy [11, 12], and others [13, 14]. In these classic treatments, p(a,t)da is used to define, at time t, the density 
of noninteracting agents with age between a and a -j- da. The total number of particles in the system at time t is thus 
n(t) = p(a,t)da. li p-{a ;n(t)) is the death rate for individuals of age a, the McKendrick-von Foerster equations 
are [11, 12] 


with p(a, t = 0) = g(a) and 


dpja.i) dp{a,t) 
dt da 


-p{a;n{t))p{a,t), 


( 1 ) 


p{a = 0,t)= / P{a;n(t))p{a,t)da ( 2 ) 

Jo 

for initial and boundary conditions, respectively. The boundary condition (Eq. 2) reflects the fact that birth gives 
rise to age-zero individuals. Note that the birth and death rates j3 and p are usually simply assumed to be functions 
of the total population n(t). 

The population dependence of /3(a; n(t)) and p(a\ n(t)) in Eqs. 1 and 2 are assumed without explicit derivation and 
it is not clear whether such simple expressions are self-consistent. Moreover, the McKendrick-von Foerster equation is 
expected to be accurate exact only when the dynamics of each individual are not correlated with those of any other. 
Therefore, a formal derivation will allow a deeper understanding of how population dependence and correlations arise 
in a fully stochastic age-structured framework. 
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Two approaches that have been used for describing stochastic populations include Master equations [15, 16] and 
evolution equations for age-dependent branching process such as the Bellman-Harris process [17-21]. Master-equation 
approaches can be used to describe population-dependent birth or death rates [11, 12, 22, 23] but implicitly assume 
exponentially distributed waiting times between events [16]. On the other hand, age-dependent models such as the 
Bellman-Harris branching process [17] allow for arbitrary distributions of times between birth/death events but they 
cannot resolve age-structure of the entirte population nor describe population-dependent dynamics that arise from 
e.g., regulation or environmental carrying capacities. 

A number of approaches attempt to incorporate ideas of stochasticity and noise into age-dependent population 
models, [3, 18, 24-29]. For example, stochasticity can be implemented by assuming a random rate of advancing to the 
next age window (by e.g., stochastic harvesting [26, 27] or a fluctuating environment [30, 31]). However, such models 
do not account for the intrinsic stochasticity of the underlying birth-death process that acts differently on individuals 
at each different age. One alternative approach might be to extend the mean-field, age-structured McKendrick-von 
Foerster theory into the stochastic domain by considering the evolution of P{n{a);t), the probability density that 
there are n individuals within age window [a, a + da] at time t [3, 32]. This approach is meaningful only if a large 
number of individuals exist in each age window, in which case a large system size van Kampen expansion within 
each age window can be applied [15]. However, such an assumption is inconsistent with the desired small-number 
stochastic description of the system. 

A mathematical theory that addresses the age-dependent problem of constrained stochastic populations would 
provide an important tool for quantitatively investigating problems in demography, bacterial growth, population 
biology, and stem cell differentiation and proliferation. In this paper, we develop a new kinetic equation that intuitively 
integrates population stochasticity, age-dependent effects (such as cell cycle), and population regulation into a unified 
theory. Our equations form a hierarchy analogous to that derived for the BBGKY (Bogoliubov-Born-Green-Kirkwood- 
Yvon) hierarchy in kinetic theory [33, 34], allowing for a fully stochastic treatment of age-dependent process undergoing 
population-dependent birth and death. 


KINETIC EQUATIONS EOR AGING POPULATIONS 

To develop a fully stochastic theory for age-structured populations that can naturally describe both age- and 
population size-dependent birth and death rates, we invoke multiple-particle distribution functions such as those used 
in kinetic theories of gases [34]. Our analysis builds on the Boltzmann kinetic theory of D. Zanette and yields a 
BBGKY-like hierarchy of equations. Here, the positions of ballistic particles will represent the ages of individuals. 

Changes in the total population require that we consider a family of multiparticle distribution functions, each with 
different dimensionality corresponding to the number of individuals. In this picture, birth and death are represented 
by transitions between the different distribution functions residing on different fixed particle-number “manifolds.” 
Processes that generate newborns (particles of age zero) manifest themselves mathematically through boundary 
conditions on higher dimensional distribution functions. 

To begin, we define 


fnixi,X2, X3,..., Xn;t)dxidx2 . ■. da;„ (3) 

as the probability that at time t, one observes n distinguishable (by virtue of their order of birth) individuals, such 
that the youngest one has age within {xi,xi -I- dxi), the second youngest has age within (a; 2 , a ;2 -I- dx 2 ), and so on. If 
the individuals are identical (except for their ages) and one does not distinguish which are in each age window, one 
can define pn{xi,X 2 ,X 3 ,..., Xn', t)dxidx 2 ■. ■ dxn as the probability that after randomly selecting individuals, the first 
one chosen has age in {xi,xi + dxi), the second has age in {x 2 ,X 2 + da; 2 ), and so on. For example, if there are three 
individuals with ordered ages xi < X 2 < X 3 , the probability of making any specific random selection, such as choosing 
the individual with age X 2 first, the one with age xi second, and the one with age X 3 third, is More generally, 

when the ages Xi_„ = x„ = (xi,X 2 ,. ■., Xn) are unordered, the associated probability density is 

Pn{^n]t) = ^fn{'^{{x^});t), ( 4 ) 

n\ 

in which T is the time-ordering permutation operator such that, for example, 7 (x 2 ,Xi,X 3 ) = (xi,X 2 ,X 3 ). Note that 
in this formulation, pn(pi.n',t) is invariant under interchange of the elements of x„. 
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To derive kinetic equations for pn{'x.n;t), we first define an ordered cumulative probability distribution 

pai pa 2 PCLn 

Qni^n: t) — / dxi / dX 2 ' ' * / ^^nfn{^n^ t), (5) 

Jo J x\ J Xn-l 

where a„ = = (oi,..., a„). Qn(s^n', t) describes the probability that there are n existing individuals at time t and 

that the youngest individual has age xi less than or equal to oi, the second youngest individual has age xi < X 2 < 02 , 
and so on. The oldest individual has age Xn-i < x„ < a„. 

We now compute the change in Qn{sLn]t) over a small time increment e: Qn(a-n + e;t + e) = Qn(^n',i) + 
where J(a„;t') = J+(a„;t') — J~{an',t') is the net probability flux at time t'. The probability 
flux which increases the cumulative probability is denoted J~^ while that which decrease the cumulative probability 
is labelled J~. Each of the include contributions from different processes that remove or add individuals. A 
schematic of our birth-death process, starting from a single parent, is depicted in Fig. lA. 

In the e: —^ 0 limit, we find the conservation equation 


dt ^ dai 

1 — 1 


J J (a^,t). 


( 6 ) 


Eq. 6 is a “weak form” integral equation for the probability density which allows us to systematically derive an 
evolution equation and the associated boundary conditions for /„(x„;t). The probability fluxes can be decomposed 
into components representing age-dependent birth and death 


J=^(a„;t) = J^(a„;t) -k J^(a„;t), 

where the birth and death that reduce probability can be expressed as 


f^ai pa2 


Jp{an]t)= / dxi / dx2- - dXnfn{^n]t)y^/3n{Xi), 

Jo Jxi J^n-l 2=1 

pCll P(^2 pOs-n ^ 

J^(a„;t)= / dxi / dX2 -- dXnfni.Xn',t)'^ Pnixi). 


I XX 


Similarly, the probability fluxes that increase probability are 


C^j + i 


n—l 


J^(a„;t) = / dxi 

Jo Jxj-i 


dXj I —if n — l (^n —1; 0 ^ ^ jJn — l ; 

Jxr,-2 


(7) 


( 8 ) 

(9) 


( 10 ) 


77- PCLl 

y(an;t)=y'/ dxi-- 

■fci../ 

ai+i pai+i 

dy / dxi+i • • 

■ / dx„/x„+i( 2 /)/„+i(xi,y,Xi+p„;t), 

(11) 

i =0 “'0 

J Xi—I J X 

i Jy 

J Xn—x 



in which x,;.j = (xj, x^+i,..., Xj), xq = 0 , a„+i = oo, and the age- and population-dependent birth and death rates for 
individual i are denoted /3n{xi) and ^„(xi), respectively. The probability flux into Qra(a„;t) arising from birth of the 
n — l individuals of age a 2 ,„ = ( 02 , 03 ,..., a„) generates an individual of age zero. Hence, a key feature of (a„; t) 
is that it does not depend on oi. 

We can now describe the fully stochastic aging process in terms of the ordered distribution function /„(x„;t) by 
using Eqs. 7-11 in Eq. 6 and applying the operator to find 


dfn{an;t) 

dt 


n p f . ^\ 71, n 

^ ^ -/-K; 0 ^ 7n(a*) + ^ 

j—1 ^ i—1 z—0 * 


l^n-\-l (y)/n+l (^2 5 Vi ^ 2 +l,n; 


( 12 ) 






4 




FIG. 1: (A) A simple age-dependent birth-death process. Each parent gives birth with an age-dependent rate j3n{a), which 
may also depend on the total population size n. Individuals can also die (open circles) at an age- and population-dependent 
rate fin{a,)- (B) Age trajectories in the upper (a > t) octant are connected to those in the lower one (a < t) through the birth 
processes. Individuals that exist at time t = 0 can be traced back and dehned by their time of birth 6;. Here, the labeling 
ordered according to increasing age. The pictured trajectories define characteristics ai{t) that can be used to solve Eq. 12. 


where ao = 0, a„+i = oo, and the total age-dependent transition rate is 

Iniai) =/3n{ai) + Hn{ai). (13) 

Note that the ai—independent source term that had contributed to the ordered cumulative (Eq. 6) does not 
contribute to the bulk equation for /„(a.„;t). Rather, it arises in the boundary condition for /„, which can be found 
by setting oi = 0 in Eq. 6. Since Q{0, 02 ,..., a„; t) = 0 and jf (a„; t) are independent of ai, the remaining terms are 


na2 rdn 

/ dx 2 --- da;„/„(a;i = 0,X2,„;t) = J^(a„;t). (14) 

Jo j Xn—1 

Further taking the derivatives of Eq. 14, we find the boundary condition 

n 

fni^l — ^2,n 7 0 — fn—1 (^2,n ; 0 ^ ^ /^n—1 (^i)- { 

i=2 

We now consider indistinguishable individuals as described by the density defined in Eq. 4. Equation 12 can then 
be expressed in terms of p„(a„;<): the probability density that if we randomly label individuals, the first one has age 
between ai and ai -|- dai, the second has age between 02 and 02 -b da 2 , and so on. The kinetic equation for can 
then be expressed in the form 


dpn dpn (a„; t) 

dt ^ doj 

j=i J 

and the boundary condition becomes 


n 

-Pn(an;t) ^ 7 n(ai) + (n + 
2=1 



00 

M«+i(2/)dn+i(a„,2/;t)d2/, 


(16) 


np„(ai,... ,af = 0,... ,a„;t) = p„_i(ai, ... ,di,... I3n-i{ai), (17) 

where the sum precludes the i = £ term and di indicates that the variable ai is omitted from the sequence of arguments 
[34]. Equation 16 and the boundary conditions of Eq. 17, along with an initial condition p„(a„;t = 0), fully define 
the stochastic age-structured birth-death process and is one of our main results. Eq. 16 is analogous to a generalized 
Boltzmann equation for n particles [34, 35]. The evolution operator corresponds to that of free ballistic motion in one 
dimension corresponding to age. However, instead of particle collisions typically studied in traditional applications of 
the Boltzmann equation, our problem couples density functions for n particles to those of n -|- 1 and n — 1 (through 
the boundary condition). 
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SOLUTIONS AND EQUATION HIERARCHIES 


Equation 16 defines a set of coupled linear integro-differential equations. We would like to find solutions for p„(a„; t) 
expressed in terms of an initial condition gn{^ — t]t = 0). However, we will see below that the presence of births 
during the time interval [0, t\ prevents a simple solution to Eq. 16 due to interference from the boundary condition in 
Eq. 17. Instead, we will obtain a solution for p„(a„;t) at time t in terms of the distribution Pnis-n — (t — to)', to) at 
an earlier time tg selected such that no births occur during the time interval (to,t]. That is, ii bi = t — ai represents 
the time of birth of the individual (see Fig. IB), we have the condition tg > biMi. The dynamics described by 
Eq. 16 are then unaffected by the boundary condition (Eq. 17) and can be solved using the characteristics at = t — bi 
indexed by individual times of birth hi. Note that any individual initially present (at time t = 0) has a projected 
negative time of birth. We can then solve pn{t — b„;t) explicitly along each characteristic and then re-express them 
in terms of a„, to obtain 

Pn (^n ,t) — Un (^n 5 tg , t)pji (tin (t to),to)“t”(n-t-I) / Un (^n ; t , t) 

Jto 

where p„+i = p„+i(a„ - (t - t'),y;t') above, and 

= U~^(aLm;tg;t')Un(aim;tg;t) (19) 


Un(^m ,t , t) — exp 


m ,t 

/ 7 «(«i - (t- s))ds 
i=l -tf 


^J■n+l(y)Pn+ldy 


dt'. 


(18) 


is the propagator for any set of m < n individuals from time t' to t. 

In the case of a pure death process where no births occur (/3„ = 0), allowing us to set tg = 0. A complete 
solution can be found through successive iteration of Eq. 18. We further simplify matters by assuming an initial 
condition that factorizes into an initial total number distribution p(n) and common initial age probability densities 

n 

g(a): /On(a„ — t;0) = p(n) J([ g(ai — t). If we further assume a death rate Pn(o) = p(a) that is independent of 

i=l 

population size, Eq. 18 can be solved, after some algebra, to yield 


n 

Pn (^n; — U ( 3.71 , 0, i) 

2=1 



p{n + k) 


t 


00 


k 


.0 s 


U{y',0-s)p,(y)dyds 


( 20 ) 


For a pure birth process where pn = 0, the second integral term in Eq. 18 disappears. In this case, we must use 
the boundary condition (Eq. 17) to successively bootstrap the solution by applying the propagator U between birth 
times. Assume a starting time t = 0 with an initial condition consisting of m individuals with corresponding ages 
a > t. The symmetry of pn(a„;t) and Un(si-n',t'',t) implies that, without loss of generality, ages can be arranged in 
decreasing order: oi > 02 > ... > am > t > Um+i > ... > a„, where the youngest was born most recently at time 
t — an > 0. If we select tg to be the moment of birth at time bn = t — an oi the most recently born (n**') individual, 
the density over all individuals is propagated forward according to 


Pn(^n,t) — Un(^mbn,t) Pn(\^n — 1 ^71,0},! an), 

where p„({a,j_i — a„, 0}; t — a„) is the initial condition immediately after the birth of the n*** individual and can be 
related to Pn-i through the boundary condition in Eq. 17. The density function thus obeys 

^ n —1 

Pn — Un{^n'i ^n‘i t)Pn—l (^n—1 0>n‘i ^ O^n) /?n—1 {Oti 

2 = 1 

Eq. 22 can then be iterated back to t = 0 to find the solution for randomly selected individuals. For the case in which 
7 „ = 7 is independent of the population size, the propagator can be separated into a product across individuals. If 
/3„ = /3 is also independent of n, the solution takes the simple form 


Pn(^in,t) gm(^m t)U(tim,d^t) U(aj^^ t) ) ld(a^ a]^) 


k-1 


k—m+l 


^=1 


(23) 
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FIG. 2: Comparison of P{a,t) (Eq. 24) derived from the McKendrick-von Foerster equation with Pm{a,t) of a fully stochastic 
pure birth process with constant /3 = 0.1. We start with = 10 individuals and analyze our quantities at time t = 10 for ages 
a < t. (A) Each of the 100 grey lines count the number of individuals younger than age a in one simulation. The solid black 
curve indicates the deterministic (McKendrick-von Foerster) solution P{a, t) = p{y, t)dy, which can also be obtained through 
P{a,t) = '^^^i'mPm{a,t). The shaded region represents the inter-quartile range of Pm{a,t). (B) Distribution constructed 
from 1000 simulations (bars) and theoretical distribution Pm{a, = 5, t = 10) (black curve). 


where bk = t — ak and gm is the initial distribution of ages for the m individuals born before t = 0. 

The above solutions for pn{aL,i]t) allow us to explicitly compare differences between the fully stochastic theory and 
the deterministic McKendrick-von Foerster model. As an example, consider the expected number of individuals at 
time t that have age between 0 and a, 


P{a,t) = f p{y,t)dy, (24) 

Jo 

where p{y,t) is found from Eqs. 1 and 2. We wish to compare this quantity with the probability Pm{a,t) that there 
are m individuals at time t with age between 0 and a. The probability Pm{n,a^t) that there are n total individuals 
of which exactly m have age between 0 and a can be constructed from our fully stochastic theory via 



The marginal probability Pm{a,t) is then found by summing over n>m'. 


OO 

^ ^ Praip^ ^ l') • 

n—m 


(26) 


The comparison can be made more explicit by considering simple cases such as an age-independent birth-only process 
with fixed birth rate /3. If the process starts with precisely N individuals, standard methods [13, 14] yields a simple 
solution of the McKendrick-von Foerster equation which when used in Eq. 24 gives P{a < t;t) = (l — e“^“). 

Substituting the pure birth solution of Eq. 23 into Eqs. 25 and 26 yields 


Pm (a, t) 


/m + N-l\ (1 

V / (1 - -b ' 


(27) 


In Eig. 2A we compare the expected value P{a, t) derived from solutions to the McKendrick-von Foerster equation 
with stochastic simulations that sample the stochastic result Pm{a,t). The fully stochastic nature of the process is 
clearly shown by the spread of the population about the expected value. Fig. 2B plots the corresponding number 
distribution Pm{5, 10). 

Finally, to connect our general kinetic theory with statistically-reduced (and deterministic) descriptions, we consider 
reduced fc—dimensional distribution functions defined by integrating p„(a„;t) over n — k age variables: 
















































7 


Pn ^ 5 0 — / ^^/c+1 ■ • ■ / dttyj (^n 5 ^) ■ (^^) 

Jo Jo 

The symmetry properties of pn{3Ln‘-,t) indicate that it is immaterial which of the n — k age variables are integrated 
out. If we integrate Eq. 16 over all ages {k = 0), and assume pn\ci = oo;t) = 0, we find 

/*oo poo 

—— =np^P{a = 0]t)-n J -fn{y)p^P {y;t)dy + {n + 1) J pn+i{y)Pnliiy^'t)^y- (29) 

Furthermore, integrating Eq. 17 over yields npn\a = 0;t) = (n — 1) /3n-i(y)Pn-iiyT^)^y- Thus, Eq. 29 can 
be written in the form 


f) r'^ 

— ={n-i)J Pri-i{y)Pn-i(.y^t)dy-nJ {Pn{y) + Pn{y))pli^Hy,t)^y + {n +1)J pri+i{y)Pnliiy^t)<iy- 

(30) 

Eq. 30 describes the evolution of the probability pP{t) that the system contains n individuals at time t and it 
contains the single-particle marginal density Pn\y', t). Upon deriving equations for Pn\y', t), one would find that they 
depend on (yi, 2/2; 0: on. Therefore, the marginal probability densities form a hierarchy of equations, as is 

typically seen in classic settings such as the kinetic theory of gases [33] and the statistical theory of turbulence [36]. 
Note that if the birth and death rates /3„ and pn are age-independent, they are constants with respect to the integral 
and Eq. 30 reduces to the familiar constant birth and death rate master equation for the simple birth-death process: 


=(n - OPn-lPn-liO - niPn + Pn)pP{t) -k (n -k l)pn+lPnll{t)^ (31) 

where p^\t) is the probability the system contains n individuals at time t, regardless of their ages. 

In general, integration of Eq. 16 over n — k > 0 age variables leaves k remaining independent age variables. The 
resulting kinetic equation for p^\ak',t) involves both (a^, j/;t) and boundary terms p^~^^\ak,ak+i = 0;t). 

These boundary terms can be eliminated by using the result obtained from integration of the boundary condition 
(Eq. 17) over n — k — 1 age variables. By exploiting the symmetry properties of the marginals we find 


dp^Pjt) ^ dp^P{ak\t) 

dt ^ dtti 

1—1 


n — k 


+ I !!— 1\+ — / Pri-i{y)Pn-Ji^k,y]t)dy 

\ ^ / ^ Jo 

^ POO 

- Pn\^k\t)^-fn{ai) - {n - k) I 7 „( 2 /)p(f+i)(afc,y;t)d?/ (32) 

Jo 


2=1 

pOO 


pOO 

-k(n-kl) / p„+i(y)p^yj(afc,2/;t)dy. 
do 


(k'] 

Each function p)i in the hierarchy not only depends on the functions in the nil subspace, but is connected to 
functions with k + \ and k — \ variables. The latter coupling arises through the boundary condition for p^P which 
involves densities p^ As with similar equations in physics, the hierarchy of equations cannot be generally solved, 
and either factorization approximations or truncation (such as moment closure) must be used. 

We now show that the fc = 1 equation explicitly leads to the classic McKendrick-von Foerster equation and its 
associated boundary condition. For fc = 1, pP^ {a]t)da is the probability that there are n individuals and that if one 
is randomly chosen, it will have age between a and a + da. Therefore, the probability that we have n individuals of 
which any one has age between a and a -I- da is npP(a; t)da. Summing over all possible population sizes n > 1 gives 
us the probability p(a, t)da that the system contains an individual with age between a and a + da: 









(33) 


p{a,t) = '^npll'>{a-,t). 

n—0 

Multiplying Eq. 32 (with k = 1) by n and summing over all positive integers n, we find after carefully cancelling like 
terms 


dp{a,t) dp{a,t) 
dt da 


OO 

=(a)pi^'> (a; t). 

n—1 


(34) 


Equation 34 generalizes the McKendrick-von Foerster model to allow for population-dependent death rates, but does 
not reduce to the simple form shown in Eq. 1. Population-dependent effects in equation for p{a,t) requires requires 
knowing the “single-particle” density function pn'^ (a; t) and subsequently all higher order distribution functions. 

A boundary condition is naturally recovered by integrating over all ages but ai in Eq. 17 and summing over all n: 


^npW(a = 0;t) = p{a 

n—1 


0,i) = 

n—2 


poo 

1 ) / Pn-1 

Jo 


(y)Pn-: 




(35) 


These equations show that the McKendrick-von Foerster equation is recovered only if both pni,ci) = y(a) and 
Pn[o) = Pia) are independent of population size. In this case, p{a) can be pulled out of the sum in Eq. 34 and 


(a;0 = P{a)p{a,t). Similarly, P{y) “ Oph\{y,t)\ P{y)p{y,t)<^y, which is the 

simple boundary condition associated with the classic McKendrick-von Foerster model. This derivation clearly shows 
that population-dependent birth and death rates cannot be readily incorporated into an age-dependent model, even 
one that is deterministic, without considering the hierarchy of population densities. 


DISCUSSION AND CONCLUSIONS 

We have developed a complete kinetic theory for age-structured birth-death processes. To stochastically describe the 
age structure of a population requires a higher dimensional probability density. The evolution of this high-dimensional 
probability density mirrors that found in the Boltzmann equation for one-dimensional, ballistic, noninteracting gas 
dynamics. However, one crucial difference is that the number of individuals can increase or decrease according to 
the age-dependent birth and death rates. Thus, the dynamics are determined by a phase-space-conserving Liouville 
operator so long as the number of individuals does not change [33]. Once an individual is born or dies, the system 
jumps to another manifold in a higher or lower dimensional phase-space, immediately after which conserved dynamics 
resume until the next birth or death event. Such variable number dynamics share similarities with the kinetic theory 
of chemically reacting gases [37]. 

Our main mathematical results are Eqs. 16 and 17. These equations show that birth-death dynamics couple densities 
associated with different numbers n and describes the process in terms of ballistically moving particles all moving 
with unit velocity in the age “direction.” The individual particles can die at rates that depend on their distance from 
their origin (birth). Particles can also give birth at rates dependent on their age. The injection of newborns at the 
origin (zero age) is described by the boundary condition (Eq. 17). 

One important advantage of our approach is that it provides a natural framework for incorporating both age- 
and population-dependent birth and death rates into a stochastic description, which has thus far not been possible 
with other approaches. In general, our kinetic equations need to be solved numerically; however, we found analytic 
expressions for pn(a„; t) when either birth or death vanishes and the other is independent of population. Furthermore, 
we define marginal density functions and develop a hierarchy of equations analogous to the BBGKY hierarchy (Eq. 32). 
These equations for the marginal densities allow one to construct any desired statistical measure of the process and 
are also part of our main results. We explicitly showed how a zeroth order equation leads to the equation for the 
marginal probability of observing n individuals in the standard age-independent birth-death processes (Eq. 31) [23]. 
The first-order equation is also used to derive a hybrid equation for the mean density p(a, t) that involves the single¬ 
particle density function pn'^ (a; t) (which ultimately depends on higher-dimensional densities {^k',t) through the 

hierarchy). Only when death is independent of population does the theory reduce to the deterministic McKendrick-von 
Foerster equation (Eq. 34) and the associated boundary condition (Eq. 35). 
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Extensions of our high-dimensional age-structured kinetic theory to more complex birth-death mechanisms such as 
sexual reproduction and renewal/branching processes can be straightforwardly investigated. The simple birth-death 
process we analyzed allows for the birth of only a single age-zero daughter from a parent at a time. We note that 
the Bellman-Harris process described via generating functions [19, 20] (which can describe age-dependent death and 
branching, but cannot be used to model population-dependent dynamics) assumes self-renewal at each branching 
event. That is, two (or more) daughters of zero age are simultaneously produced from a parent. Such differences in 
the underlying birth process can lead to qualitative differences in important statistical measures beyond mean-held, 
such as hrst passage times [21]. The branching/renewal process, as well as sexual reproduction, requires nontrivial 
extensions of our kinetic theory and will be explored in a future investigation. 
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